AD-A036  391 

UNCLASSIFIED 

I °F'l 

^AOTSMI 


STANFORD  UNIV  CALIF  INFORMATION  SYSTEMS  LAB  F/«  12/2 

New  RESULTS  IN  2-0  SYSTEMS  THEORY#  2-D  STATE-SPACE  MODELS  - REA— ETC (U) 
1976  B LEVY#  S Y KUN6#  M MORF  F44620-69-C-0101 

AFOSR-TR-77-0370  NL 


DATE 
■ FILMED 

5-77 


elements.  The  l-1)  realization  problems  have  been 
veil  studied  and,  given  any  transfer  function.  It 
Is  well  known  that  the  realization  can  be  readily 
found  in  certain  standard  (e.g.  controller  canoni- 
cal) forms  [11].  For  the  realization  of  a 2-D 
transfer  function,  a major  difference  is  that  two 
types  of  dynamic  elements  are  needed  - "horizontal 
delay  element"  (z-3)  and  "vertical  delay  element" 
(a:-3).  Now  an  important  problem  is  that  of  how  to 
use  2-D  dynamic  elements,  multipliers  and  adders  to 
realize  a 2-D  digital  filter  with  the  transfer 
funecion: 
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Figure  3.  2:  2-D  Controller  Form  Realization 
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We  can  do  this  in  two  steps.  First  we  rewrite 
(3.1)  in  a rational-gain  representation,  i.e. 


H(z~\(i)_1) 


.■I  b.((i)_1)  z*'1 

1=0  1 

[ a (of1)  z_1 
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(3.2) 
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Without  loss  of  generality,  we  can  assume  a 
and  we  denote 

a0(u_1)  A i +T0((j‘1) 

Thus,  using  the  1-D  realization  technique,  we  write 


The  realization  is  almost  achieved:  in  addi- 

tion to  the  n horizontal  delay  elements  we  need 
only  n vertical  delay  elements  to  implement  the 
feedback  gains  {a{(w-3),  i=0,  l...m}  and  m other 
vertical  delay  elements  to  implement  the  readout 
gains  {bi(oj'1),  i=0 , l...m}.  Thus,  the  complete 
realization  shown  in  Figure  3.2  requires  only  n+2n 
dynamic  elements. 

This  realization  is  a standard  (canonical) 
one;  its  structure  is  very  simple  and  it  involves 
only  real  gains.  Note  also  that  we  need  fewer 
dynamic  clcnents  than  the  implementations  in  [9], 


State  Space  Model  Representation 

As  remarked  in  Section  2,  circuit  Implemen- 
tations with  delay  elements  z-3-  and  u'3  are  in  a 
one-to-one  correspondence  with  state  space  models 
of  Roesser's  type.  The  output  of  the  z~3  delays 
are  the  horizontal  states  and  the  outputs  of  the 
w-3-  delays  are  the  vertical  states. 

Thus,  the  implementation  of  Figure  3.2  can  be 
transformed  readily  into  the  following  state  space 
model : 
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state  space  models  can  be  tound  In  l XU  J . 

Hardware  Design  of  2-D  Digital  Filter 

The  idea  of  using  tuo  types  of  dynamic  ele- 
ments is  not  very  abstract;  it  is  very  natural  in 
delav-dif for on dal  systems.  However,  before 
considering  its  practical  application  to  image  sys- 
tems, two  remarks  have  to  be  made: 

1)  Because  the  "spatial"  dynamic  elements  seem 
unimplemer.cable,  (except  as  index  operators  in  a 
digital  computer,  for  example,)  we  can  replace  them 
by  time-delay  elements. 

2)  In  order  to  have  a f inite  order  descript  ion , 
we  shall  only  consider  a bounded  frame  system,  i.e. 
we  assume  that  the  picture  frame  of  Interest  is  an 
MxN  frame  (with  vertical  width  M and  horizontal 
len.gch  N)  . 

Note  that  in  order  to  use  time  delay  elements 
ve  need  first  to  find  a way  to  code  a 2-D  spatial 
system  into  a 1-D  (discrete-time)  system  and  vice 
versa.  Thus  we  shall  propose  the  following  imple- 
mentation of  a 2-D  filter; 

i)  The  input  scan  generator  codes  the  2-D 
spatial  input  into  1-D  (time)  data  according  to 
the  mapping  function  t(.,.) 

- 1H  + JN  (3.4) 

where  M and  N are  relatively  prime  integers. 

it)  A 1-D  (discrete-tine)  digital  filter 
processes  the  1-D  data  generated  by  (i) . This 
subsystem  is  implemented  by  replacing  z-i  by  6 , 
ta”l  by  A in  a 2-D  circuit  realization  (e.g.  2-D 
controller  form) . 6 and  A are  chosen  as 

M 

6 *=  D = M-units  delay  element, 

N f3'5) 

A » D “ N-units  delay  element. 

iii)  The  output  frame  generator  decodes  the 
1-D  (discrete-time)  output  of  the  1-D  digital 
filter  described  above  into  a 2-D  (discrete-spatial) 
picture  according  to  the  inverse  mapping  of  (3.4). 

(i(t) , j(t))  = (Ptmodtt,  [t-  (PtmodM)N]/M)(3.6) 
where  P is  the  unique  integer  such  that 

PN  - QM  =>  1 and  0<P<M  (3.7) 

Verification:  Let  us  note  the  1-D  (discrete- 

time) output  will  be 

i/(D)  - H(D)  U(D) 

- H(z  3,w  ^)u(z  *) 

-1  _M  -1  N 
z =D  ,w  -D 

r r — i -j  | 

- Z Z y . , z 0)  J 

1 J Iz-WVW*  (3.8) 

where  (ytj)  represents  the  2-D  (discrete-spatial) 
output  data  field.  Note  also  that 

i/(D)  A Z t f D-t  t (3.9) 

t c 


yt-  l y,  . 

(i,j):iM+jN  - t 


(3.10) 


Since  the  system  is  a causal  system, 

y±i  "0  if  i,j  < 0 . (3.11) 

Let  ur  consider  only  the  integer  t with 
t = iM+JN  , i < N , j < M 
then  (3.10)  and  (3.11)  give 


since,  for  this  special  case,  the  summation  set  of 
(3.10)  contains  only  one  nonzero  point.  There- 
fore, we  will  obtain  a bona  fide  output  picture 
inside  the  MxN  frame. 

This  2-D  image  scanning,  and  display  system  is 
not  as  complicated  as  it  looks,  it  can  be  simple: 

Example;  Problem:  Design  a 2-D  digital  filter  for 


H(z_1 ,u-1) 


1+0.  3z  X+0 . 2(J_1+0.  lz_^U)_^ 


for  a frame:  MxN  = 100x101.  Assume  D=*0.0l  ms. 
Solution.  (i)  ISO 

In  this  special  frame  (with  ti=H+l),  the  input 
scanning  generator  is  indeed  very  simple,  as  shown 
in  Figure  3.3. 


lms/iint 

Scanning  time:  0.01  ms/pixel  = 1 ms/line 

Scanning  angle:  45° 

Fig.  3.4  Input  Scan  Generator  & Output  Frame 
Generator 

(ii)  1-D  digital  filter 

Constructing  the  2-D  realization  of  Figure 
3.2  and  then  replacing  z“l  by  6 and  by  A we 
have  the  1-D  realization  shown  in  Figure  3.4 


1.01  ms  delay  element 
1.00  ms  delay  element 
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The  output  frame  generator  does  the  reverse  of 
the  ISO,  displaying  the  picture  instead  of  scanning. 

Dimensionality  of  Clobal  State 

Considering  a bounded  frame  (MxN)  system,  it 
is  interesting  to  knew  the  dimension  of  the  global 
state  (or  initial  conditions)  needed  to  process  the 
MxN  "future"  data  field.  Since  vertical  states 
convey  information  vertically,  all  the  vertical 
states  along  the  X-axis  are  necessary  initial  con- 
ditions and  their  dimension  is  nil.  Similarly,  all 
Che  horizontal  states  along  the  Y-axis  are  neces- 
sary initial  conditions  (with  dimension  nM)  since 
they  convey  informaciou  horizontally.  Therefore, 
in  the  bounded  frame  case  a total  number  of  mN+nM 
are  needed  Co  summarize  the  "past"  information. 

This  very  same  idea  can  be  used  again  from  a 
computational  point  of  view.  Indeed,  the  number  of 
required  storage  elements  for  recursive  computa- 
tions is  also  equal  to  nM+mN  if  initial  conditions 
are  not  zero.  However,  the  Initial  conditions  are 
often  zero,  then  chesize  of  storage  required  can 
be  reduced  to  nH  (resp.  nM)  by  storing  the  updated 
data  row  by  row  (resp.  column  by  column).  No 
storage  is  needed  for  the  rest  of  the  initial  con- 
ditions - nM  horizontal  states  (resp.  ndi  vertical 
stares)  - since  they  are  assumed  to  be  zero.  This 
Is  consistent  with  the  result  of  Read  [12]  derived 
from  a direct  polynomial  approach. 

Another  interesting  observation  concerns  the 
dimension  of  the  1-D  digital  filter  contained  in 
our  2-D  digital  filter  design  discussed  above. 

Since  it  needs  r.  M-unit-delays  and  n N-unit-delays , 
the  corresponding  1-D  state-space  has  also  a dimen- 
sion equal  to  nM+mN.  Note  that,  despite  the  high 
dimension  of  the  corresponding  1-D  filter,  its 
high  sparsity  is  very  encouraging  for  further 
studies. 

In  short,  our  studies  on  the  dimensionality  of 
2-D  global  states  have  reached  a consistent  conclu- 
sion from  either  theoretical  or  practical  approaches. 


4.  Global  and  Local  Controllability  and 
Observability 

For  reasons  of  space  we  deferred  this  Section 
to  [18],  part  II. 

5,  Modal  Controllability  (Observability)  and 
Minimality 

In  the  1-D  case,  the  relative  priaeness  con- 
cepts could  also  be  used  to  define  controllability 
and  observability.  In  [16]  Rosenbrock  proved  that 
A,  B was  controllable  if  and  only  if  zI-A,3 
were  left  coprime. 

C,  A was  observable  if  and  only  if  C,  zI-A 
were  right  coprime. 

This  approach  can  be  generalized  very  easily 
to  2-D  systems  and  will  also  provide  a definition 
of  ainlmallty. 

Definition  5.1  Let  H(z,w)  » VT_1U  where'  V,  T,  U 
are  2-D  plyr.oatal  matrices.  It  is  a minimal 
description  of  H(z,u)  if  and  only  if 

( 


V/'f  are  rTghl "coprIme”nnfTr,  U JWPUtft  co 

prime . 

This  amounts  to  requiring  that  there  is  no 
cancellation  in  the  2-D  transfer  function  H(z,ei).Ln 
[18]  part  I we  also  provide  the  important  property' 
that  if  (V,T,U)  and  (V].,Ti,Ui)  are  two  minimal 
descriptions  of  H,  |t|  = [TjJ.  We  also  presented 
an  algorithm  to  extract  the  greatest  common,  right 
(left)  divisor  of  two  polynomial  matrices,  which 
enables  us  to  find  a minimal  description  of  H from 
a nonmlnlmal  one. 

Define  A(z,u)  A [ ( ^ ) - A]  = A^^z.u)  . 

Then,  in  the  state  space  description  case 
H = CA(z,(i))--1b  is  minimal  if  and  only  if 

(A  (z,u)  , B]  are  left  copriz:e  (5.1) 

and 

[C ,A(z,u) ] are  right  coprime  (5.2) 

Definition  5.2  (i)  A,B  is  modally  controllable 

if  (5.1)  holds. 

(ii)  C, A is  modally  observable  if 
(5.2)  holds. 

These  definitions  are  clearly  connected 

to  minimality  but  the  state  space  signifi- 
cance of  controllability  and  observability  disap- 
pears in  this  formulation.  Tills  is  why  we  shall 
give  now  an  equivalent  state-space  characterization 
of  the  notions  of  modal  controllability  and  obser- 
vability . Another  consequence  is  that  for  a 
single  input-single  output  system,  if  H(a,a)  » 

k and  if  b and  a are  coprime  with  3 a = n 
a(z,ai)  , z 

3ua  = m,  then  if  CApq(z,w)  B is  a minimal  realiza- 
tion of  H(z,u)  we  must  have  |A(z,uj)|=  a(z,!c)  and 
hence  p = n and  q = m. 

Hence  the  validity  of  our  definition  of  mini- 
mality of  a state-space  model  will  depend  on  our 
ability  to  realize  a transfer  function  of  order 
(n,n)  with  n+m  states.  This  problem  was  considered 
in  Section  6 of  [18],  part  II. 

A consequence  of  the  relative  prineness 
criterion  for  2-D  polynomial  matrices  given  in 
[18],  part  I is  that  C and  A(z,w)  are  right 
coprime  if  and  only  if 

tank  p\(z,w)j  = n+ri 

for  any  generic  point  (£j,  C2)  any  irreducible 
algebraic  curve  appearing  in  the  decomposition 
of  V,  the  algebraic  curve  defined  by  |A(z,u)|  = 0. 

It  is  to  be  noted  that  the  rank  is  considered  over 
the  field  K(Ci,  g 2)-  A proof  of  this  is  given  in 
[18],  part  II,  along  with  some  illustrating 
examples. 


6.  Minimality  of  State-Space  Model 

It  is  shown  in  the  last  section  and  in  [18], 
part  II,  that  only  a state-space  realization  with 
order  (n,m)  - i.e.  the  same  order  as  the  transfer 
function  - can  be  both  modally  controllable  and 
modally  observable.  Now  the  question  is  whether 
such  a realization  exists  at  all. 

The  best  way  to  prove  the  existence  of  such 
realization  is  by  construction.  Note  that,  in  the 


s 

J ' 


2-D  state-space  model,  the  particular  transform 
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enables  us  to  change  the  basis  of  the  state-space. 
The  matrices  {a,B,C,D}  are  transformed  to 


T A T"1 
CT-1 


B - TB 
D-D 


In  fact,  it  is  more  convenient  to  work  with  a ca- 
nonical form  under  the  "similarity  transform" 
defined  by  (6.2). 

In  the  1-D  case,  all  minimal  state-space  model 
can  be  transformed  to  the  controller  canonical  form. 
Similarly,  almost  all  [18]  2-D  state-space  model 
can  be  transformed  to  the  following  modal  controller 
form  [A,  B,  C,}  (assuming  D = 0) 


_A2i  : V-i-o 


B = |e'  ] 

c - 


A(z,w) 

b" 

et 

-C 

0_ 

where  Z . , a.  ,b.  ,e.  were  defined  in  (3.3)  and  the 
entries  of  A^2  an<^  A21  are  to  be  chosen  such  that 

det [A(z,u) ] = a(z,w)  (6.4) 


b(z,ta) 


It  is  easy  to  check  that,  in  (6.4),  the  co- 
efficients {a^,0<i<n}  and  {aoj  ,0_<J_<m}  have  already 
been  matched.  Similarly,  in  (6.4),  the  coeffici- 
ents (biQ.  0£l£n)  and  {bgj  , 0£j<m}  have  already 
been  matched.  Therefore,  only  2nm  coefficeints 
(a^j  j l<i<n  , l£j£n}  and  [bij  , l£i<n  , l£j£m}  are 
to  be  matched.  In  other  words,  there  are  totally 
2r.m  (nonlinear)  equations  to  be  satisfied.  Coinci- 
dently,  the  number  of  free  parameters  in  matrices 
A12  and  A21  is  also  2nm.  Therefore  it  is  natural 
to  conjecture  that  a solution  (or,  more  precisely, 
a finite  number  of  solutions)  should  always  exist. 

Now  let  us  examine  the  plausibility  of  this 
conjecture  by  taking  a low-order  example. 

Example  6.1  (1,1)  order  case 

For  ease  of  notation,  let  A^  ” a » A2].  = B. 
Also  (without  loss  of  generality)  let  us  assume 
that  b^o  r 0 (otherwise,  we  may  have  to  use  another 
canonical  forn) . Then  (4)  becomes 

zoH-a-,.  z+a.n- u-d3  “ zuthagi- z+a^g*  w+all  or  equiva- 
Q . lently 

03  “ '*11  (6.6) 

and  (6.5)  becomes 

b01Z+b10^(a10b01+a01b10a+b013)  • b01Z+b10“*bll 


b01B+b10a 


bll"a01b10"a10b01 


a “ 2b10  (b10-a01b10-a10b01±/(bira10b10-310b0l‘ 

““  -4allb01b10r 

B- 

p a . (6.8) 

Therefore,  the  existence  of  (1,1)  order  state- 
space  model  has  been  proved  by  construction . 

Unfortunately,  (6.4)  and  (6.5)  usually  give 
a set  of  2nm  nonlinear  equations;  therefore  the 
solution  may  not  always  be  in  real  numbers.  For 
realization  with  real-gain  constraints,  we  often 
need  a realization  order  higher  than  n+m.  To  show 
that  an  (n,m)  order  real-gain  realization  may  not 
exist,  it  is  easiest  to  work  on  an  example. 

Example  6.2  The  problem  is  to  show  that  there  is 
no  (1,1)  order  real-gain  realization  for  the 
transfer  function  , 

z + to 

zoo  - 1 * • 

Solution:  Let  us  assume 


~0  cf]  P e 

, B - 

_fi  Oj  [_f_ 


c - [g,  h] 

Since  an  = -1,  B “ ct 
and  (6.5)  becomes 

fhz  + egw  - (eha 

or  equivalently, 


(6.9) 

Then  (6.4)  is  satisfied, 


eha  * + gfa 


+ gfa) 

fh  - 1 
eg  = 1 
fa  = 0 


(6.10) 

(6.11) 

(6.12) 

(6.13) 


Now,  (6.13)  x hg 
gives 


(6.11)  a g a - (6.12)  x hV 

2 2-1 
ga+ha  =0. 


(6.14) 

Since  (6.14)  has  no  real  number  solution,  no  (1,1) 
order  real-gain  realization  exists,  (e.g.  f = h = 
e **  g 1 , a = -B  = Aly)  . O 

In  the  practical  aspect,  real-gain  realiza- 
tions are  much  more  desirable  than  complex 
realizations  because  the  former  are  much  easier  to 
physically  implement.  Therefore  our  (2m+n)  order 
real-gain  realization  (cf.  Section  3)  are  justi- 
fied to  be  practical  and  low  .order  realizations. 
Indeed,  for  the  transfer  function  in  Example  6.2, 
the  minimal  real-gain. realization  {A,  B,  C}  can  be 
obtained  by  our  realization  method; 

f0  | -1  0l 


[1 ! 01] 
[1|01] 


(6.15) 


Since  b{Q  -f  0.,  (6.6)  and  (6.7)  have  solutions 
*E.  Sontag  (Iniv.  of  Florida)  independently  arrived  at 


Special  Transfer  Functions 

In  designing  digital  filter,  the  transfer 
function  may  be  intendedly  chosen  in  a certain 
form  for  the  purpose  of  an  easier  and/or  better 


the  same  conjecture  recently  (private  communication). 


realization.  Therefore,  It  Is  worth  mentioning 
that  some  special  types  of  transfer  function  can  be 
easily  realized  in  (n+m)  order  real-gain  realiza- 
tions. There  are  two  important  special  types  of 
transfer  functions: 

. 1.  with  separable  denominator, 

il.  with  separable  numerator. 


Let  us  first  consider  the  separable  denominator 
case.  Assuming 

-1  -1  -1  -1  ' <6'16> 
-1  -1,  b(z  ) b(z  ,cj  ) 

h(z  ,u  ) - -rf-rr;  - -v  -r-r-  . 0 

a(z  ,o)  ) (z  ,B(u>  ) 0 0 

n n 


l l V z"1  m'i 
i-0  j-0 


'ij 


(Oq-kIj^z  1+. ••+C1nz  1+...+0bu  m) 

then  its  circuit  realization  is  shown  in  Figure  6.1 


Secondly,  let  us  consider  the  separable 
numerator  case,  which  is  to  say  a system  with 
transfer  function. 


H(z  \«  *)  - lH(z  1,u  1)]  1 


ct(z~l)  6 (a)”1) 


l l b 

i-0  j-0 


ij 


-hri 


Z b) 


(6.17) 


At  first  sight,  it  seems  quite  difficult. 
However,  in  actuality,  the  realization  can  be 
readily  obtained  bv  using  the  inversion  rule  by 
Rung  [17].  .More  precisely,  to  realize  the  inverse 
system  of  Figure  6.1,  we  first  note  that  the  path 
"input  — — input"  is  a "feed  through" 

path  (i.e.  a path  connecting  input  and  output  with 
only  constant  gains).  The  second  step  is  to  Invert 
all  the  gain3  and  reverse  all  the  arrows  on  the 
path  (in  our  case,  replace  by  bgg“l) . Lastly, 
change  signs  of  the  gains  of  the  branches  which  are 
entering  this  path.  These  steps  complete  the 
realization  of  H"l(z”i,  or*),  In  [IS],  part  II  the 
implementation  is  given. 


Remark:  In  many  design  problems  the  constraints 

on  numerator  are  much  weaker  than  on  denominator, 
hence  this  second  form  seems  to  have  higher  poten- 
tial in  practical  applications. 
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Anr>endtx:  ?-n  Levinson  Algorithms:  The  following 

set  o:  results  were  motivated  by  the  problem  of 
determining  stability  of  2-D  recursive  filters. 

In  the  1-D  case  the  connections  between  stability 
orthogonal  polynomials  and  the  Levinson  recursions 
are  by  nou  well  known. 

In  the  2-D  case.  Shank  conjectured  that  the 
least- squares  inverse,  say  b(z,i±<),  of  an  unstable 
2-D  polynomial  a(z,w)  (of  degree  n in  z,  m 
in  u)  is  stable,  i.e., 

n ,m  . i 

t b z or  =■  b(z,u)  ~ l/a(z,u)  J (A.l) 


i,j  = 0,0 

where  b(z,u)  minimizes 


a(z,u)b(z,ui)|| ' 


e^Q  A [1,0,. ...0],  [1  x (2n  + 1)  (m  + 1)] 

- = [b00’“-,b0m’b10*''"bnm] 

and  the  Toeplitz  block  Toeplitz  matrix  CL  contain- 
ing the  coefficients  of  a(z,o)  such  that  . [Cl  b] 
is  the  vector  of  the  coefficients  of  the  product 
polynomial  [ a(z,u)b(z,u)  ] . Then  1j  is  given  by 
the  solution  of 

« b = [a'alb  = ^ aQ0  (A. 3) 

By  applying  the  Levinson  (LR'.J)  recursions  [19],  for 
block  matrices  developed  by  Robinson  and  Wiggins 
[ T03>  to  (A.  3)  b can  be  obtained  from  the  first 
column  of  the  block  solution  9,  15  = [ Im,0, . . . ,0] 1 . 
Or  with 

U A [ l,(i), . . . ,um] ' 


b(z,u)  - Z b.(z)u  = u'b(z)  (A. A) 

i=0  x 

- [u'  ,(o'z,...,ii>'zn]B  £l  = u'BCz)^  . 

Using  the  property  of  the  LRW  recursions  that 
|3(  z)J  has  its  roots  inside  the  unit  circle,  the 
g.c.d.  of  (bj_(:)},  which  divides  |B(z)|  contains 
a subset  of  these  roots. 

We  can  therefore  conclude  that  the  nonprinitlvc 
f ac  tors  of  b < z --the  contents  in  z and  u--are 

indeed  stable. 

However  Genin  and  Kanp  [21]  proved  that 
b(z,tz),  therefore  the  primitive  factors,  are  in 
general  not  stable  for  n,m  > 1. 

A 2-D  I.rvinsop  Algorithm:  Genin  and  Kamp  developed 

a 2-D  generalization  of  the  orthogonal  polynomials 
on  the  unit  circle.  We  give  here  an  equivalent 
recursion  in  the  time-domain  using  a stochastic 
framework  (sec,  e.g.,  [19]). I 

We  consider  a finite  window  of  a scalar  2-D 
stationary  stochastic  process  (yjj,  i e[0,n],j  €[0,m]] 
with  zero  mean  and  covariance 

E(yl |y_hk  \ ” _r_l-h  ,J_-k  ' 

T AT”£aggero"'r"  (MH  ) independently  also  developed 
such  recursions  (private  communication). 


pc fine  (n,m) 

*J  - Cyj,.n"”’yj,0] 

as  the  column  of  the  data  array,  0 < J < n 

and 

,/n,m)  ^ [v(n,n)i  ]t  f 

the  data  array  scanned  column  by  column.  Now  the 
covariance  of  ny  (n > ni)  is  given  by 

EjyO*.*),, («,«)')  _a<n.o)  (A>3) 

where  R(n>n)  is  a (n+l)  by  (n+1)  block  Toeplitz 
matrix  with  Toeplitz  block  entries  R;_^  = [a^>~)i 
of  size  (mi-1)  by  (mfl)  and  R_jc  =>  R^, 

(RkUj  “ rk,j-L-  Now,  let 

A 

y(n,t;n,m)  = E(y(n,i)jy(r,s)  : 

(0,0)  < (r,s)  < (n,m)  , (r , s)  i (n,f.)), 

0 < l < m (A.  6) 

^ ^ (0,ra)  (n,n) 


(0,0) 


y(n,t;n,tn)  = - y(n’m)  'a^*^  , (A. 7) 

where 

a(n!$  - [a£n’™)'(n>U>--->ann’n)('',U]  1 

and 

t im-t+l  position 

Qq  (a, (,)  — [x, ...  ,x,0,x, . ..  ,x] 

So  that  if 

V £ 9 c[y(n,n;n,m) ,...,y(r,s) ;n,n) J 

ex.  ex 

y(n,t;n,m)  ■=  y(n,^.;n,m)  - y (n,-l;n,m) 


~(n,m)’  _ _ v(n,m)  *a(n,m) 


a(n*n)  = e.©i  1 + [a<n’m>:  ...!a<n’^] 

c -lw  mfl  L (n,tn)i  J (n.f,)!  i (n,0)J 

(A. 9) 

Note  that  diagonal  entries  of  top  block  of 
Q(n,m)  equal  unity,  and  © denotes  Kronecker 
product.  Similarly  we  can  define 

y (k,m;n,m)  = E(y (k,m) | y (r , s) : 

(0,0)  < (r,s)  < (n,m) , (r,s)  ? (k,n) } = (A. 6') 

, . , 0 < k < n (A. 7’) 


X.  ^ I X.  X.  X. 

’ =*[y (n,m;n,m)  , ...  ,y  (k,m;n,m)  ....  ,y(0,m;n ,m)  ] , 

y(k,tn;n,m)  - y(k,m;n,m)  - y(k,m;n,m) 


- . v(n.®),a(n.»)  (A.8') 

r r 

where 

r tH-l®-l  la(n,ra)|  "*a(k,m)  |”*iG(0,nO  VA*9  ' 

Note  that  diagonal  entries  of  block  composed  of 


ith  row  of  ictl  block  entry  of  <2; 


equal  unity. 


Also,  the  first  columns  of  U^n,,n)  and  of  (j(n*m) 
are  the  same.  Then  by  (A. 5),  (A. 8)  and  (A,8'J 

E(y(niU)[v^‘ay>  5<n  ’ 1 m)'  ] }=«(n  ’ 1 n)  [ a <n ,ra)  ,a^n,m)  ] 

and  also,  by  (A. 6)  and  (A. 6') 

E[t/^n ,n)[?^n,ln> ' ,v(n,m) ] ) = [£(n»m)>£(n>m)] 

£^1’^Dias(e(n-i)  .njn.E)  ) 0 [ (mf  1)  x (rri-1)]  ^ 

e(k,m;n,m)  > 0 
€(n,-t,;n,m)  > 0 . 

£(n>m)  a ^ 0Diag(e(n,m-i;n,m)  } , 

hence 

a*n,n)[a^n,t!^,a^n,m^]  = [e^n,m),e^n,!n)]  (a.io) 

are  the  2-D  Levinson  ecuatlons . therefore  we  have 
n+mfl  auxiliary  solutions.  Also,  the  first  column 
of  C£n,u)  (or  of  Q(n>ra)  since  they  are  the  same) 
corresponds  to  the  2-D  causal  estimate  of  y(n,m) 
given  y (i, j)  : (0,0)  < (i,j)  < (n,m),  i.e.,  it  is 
the  ++  predictor  of  y(n,m)  (one  quadrant- 
predictor).  The  last  column  of  a£n>m)  gives  the 
— f predictor  and  the  last  column  of  gives 

the  + - predictor. 

The  Levinson  Recursions:  First  define 

[ J ],  . A /„  if  i=m'j  , 0 < l,j  < m 
1 mJij  - \0  else  ’ - ,J  - 

Now,  observe  that  the  following  reorderings  hold: 

On®' ’ <A*U> 

hence 

I(rvt-l)  (mfl) 

- CJn®V^m)-£cn,n,)J  * 

so  that 

m(a,a)(j  )[a(n’n0,a(n,n')]  [e(n’m),e(n,m)] 

n — ' m r c r c 

and  we  multiply 


£*(n,m)  . (J  gj  )£<n»m>j 

v m r n ’ 

e*(n,m)  - (j  ©j  )e(n,0)j  . 

c n ^ nr  c m 


a(n.»)[a*(n,m)a*(n,m)]  . [£*(n,^(n,»)](A<12) 


*(n,m) 

£r  “ Jn  Dia8(e(o-i,n;n,m)  )Jn©e_rf_1  = 

=■  E*(n'm)©e  , 

(we  have  used  (A  ® B)  (C  ® D)  =*  AC  © BD)  . Similarly 
£ ^n,m)  = Tirt-l©Jm  DiaSfe(n,m-i;n,m)  Jj^  = 


*(n,m) 


^nri-l^-c 


Now  define 


«f’”)  a ,...,r  la'"'"0  (a. 13) 

) '"l  <a-14> 

Now,  the  2-D  Levinson  recursions  can  be  described 
as  follows.  Increase  in  n:  n -»  n+1,  m m. 


— 

°mf  1 

^(n,n)' 

C 

c 

a*(n,m) 

m 

e(a,m) 

c 

£*(n,m) 

c 

_°m+l 

^(n,™) 

- c 

(/(n*m)  = J C*n’m)j  ) 


me  m 


Now,  let 


— (n+l ,m) 
uc 


*(n,m)-l,(n,m)|  (n,tn)-l 


(A. 15) 


A<n,n)  _ E(n,m)  _ ^(n,in)E*(n,in)-ltf(n,n>) 


on  the  right  and  denote 
*(n,n)  _ .. 


(j  ©j  )a(n,m)j  , 

N n ^ nr  r n 


a*(n,m)  . (J  gj  )C(n,m)J 
c ' m7  c it 


a(n+l,n)5(n+l,m)  @J 


1 ^ mfl 


(A.  16) 


and  the  diagonal  of  the  top  block  of 

equals  Dlag(e(n+l,m-i;n+l,m)  ) * so  that 
is  just  obtained  by  re-normaliting  the  columns  of 


— (rH-1  ,m) 
c * -1 

Note  1:  e (n+l,k;tnfl,n)  i 0 otherwise  deleting 

the  column  and  kch  row  of  t^P+li®)  , we 

could  find  cjc : R,.  + ,ra'?jc  *=  0.  But 

is  a covariance  and  this  would  mean  that  the 
estimation  problem  is  singular. 

Note  2:  Similarly  e (n+1 ,k;n+l ,m)  J 0 otherwise 

there  would  be  c:  jft(n+l,ra)c  = 0.  Also 


(n+1 ,k;n+l,m) 
(n+l.ra).  = o 


(r+l,n) 


and  let 


. (r+1 ,m) 


_(n+l,m) 

G_ 


Im4-1.. 

£(n.m) 

o: 

_ r 

)a(n,m) 

Pr 

(A.  17) 

. « 2 ■ 3 E 

H-  O ?L 

■j  c*  3 O **■ 

/pig 

B 3 “ ^ «. 

iii  v;  >— 

**  V-  a 8 

o b « s >?■ 


(tvfl.m) 


ri>“'  - |a(f1'm)  I avlrri}mM  , (n+2  columns) 

L | J (A.WJ 

re  is  the  first  column  of  G^n+*,m^. 

ci  c 

To  obtain  the  recursions  for  an  increase  in  m, 

just  have  to  reorder  in  blocks  of  size 

l)x(n-rl)  then  the  roles  of  m and  n are 


-(irt-l,ra)J 


we  just  have  to  reorder  R'n>m'  in  blocks  of  size 
(n+l)x(n-rl)  then  the  roles  of  m and  n are 
exchanged  as  well  as  (2C  and  CL r and  we  can  use 
the  same  recursion  as  the  one  just  described. 

This  version  of  the  recursion  enables  us  to 
increase  n and  m separately,  instead  of  the 
scheme  proposed  by  Genin  and  Kamp  where  m = n. 

Note  3:  The  inversions  required  by  these  recursions 

have  additional  structure,  i.e.,  the  matrices  are 
typically  non-Toeplitz , but  sums  of  products  of 
Toepl'itz  matrices.  One  can  take  advantage  of  such 
a structure  by  using  generalized  Levinson  recursions 
[22]  to  find  a representation  of  the  inverse  of  such 
natrices  also  in  terms  of  sums  of  products  of  Toeplitz 
matrices.  Expressions  with  Toeplitz  matrices, 
since  they  are  related  to  convolutions,  can  be 
evaluated  using  Fast  Fourier  Transforms  (FFT' s) . 
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